options(jupyter.plot_mimetypes='image/png')
options(repr.plot.width=8, repr.plot.height=6)

Conceptual

1. 책의 식 (5.6)을 유도할 것.

이건 하나하나 치는 게 너무 짜증나서 사진으로 대체하도록 하겠습니다.

library("IRdisplay")
display_png(file="picture.png")

2. (a)

bootstrap의 첫 번째 자리에 원래 모집단의 j 번째 값이 발견될 확률은 $\frac{1}{n}$ 이다. 그러므로 발견되지 않을 확률은 $1-\frac{1}{n}$ 이다.

2. (b)

bootstrap은 복원추출을 허용하므로, 역시 동일하게 $1-\frac{1}{n}$ 이다.

2. (c)

bootstrap sample의 각 자리에서 원래 모집단의 j 번째 값이 발견되는 사건은 모두 독립이다. 그러므로 각 사건이 일어나지 않을 확률을 모두 곱하면 어느 자리에서도 원래 모집단의 j 번째 값이 발견되지 않을 확률이 된다.

$(1-\frac{1}{n})^n$

2. (d)

n = 5 일 때, 원래 모집단의 j 번째 값이 bootstrap sample 안에 있을 확률은, 전체 확률인 1에서 하나도 관찰되지 않을 확률을 빼서 구할 수 있다.

$1-(1-\frac{1}{5})^5 = 1-(\frac{4}{5})^5 = 67.2\%$

2. (e)

n = 100 일 때,

$1-(1-\frac{1}{100})^{100} = 1-(\frac{99}{100})^{100} = 63.4\%$

2. (f)

n = 10000 일 때,

$1-(1-\frac{1}{10000})^{10000} = 63.2\%$

2. (g)

plot을 아래에 그립니다.

x = 1:100000
pr = function(x){
    return (1-(1-1/x)^x)
}
plot(x, pr(x))

원래 모집단의 j 번째 값이 발견될 확률이 63.4% 근처로 빠르게 수렴하는 것을 알 수 있습니다.

2. (h)

책의 코드를 보고 코멘트를 해 보겠습니다.

# 우선 store 라는 이름의 벡터를 만듭니다. 이 벡터는 10000개의 빈 공간을 가집니다.
store = rep(NA, 10000)

# store의 첫 번째 값부터 마지막 값까지 채워 넣을 겁니다.
for (i in 1:10000){
    store[i] = sum(sample(1:100, rep=TRUE) ==4) > 0                
    # sample(1:100, rep=TRUE)
    # 1부터 100까지의 값들 중에 복원추출을 통해 무작위로 뽑은 
    # 100개의 값들로 이뤄진 벡터를 구합니다

    # sample(1:100, rep=TRUE)==4
    # 그러고 난 다음에 벡터의 각 값들이 4와 같은지 비교합니다.
    # 이렇게 하면 원래 숫자였던 벡터의 값들이 TRUE 혹은 FALSE가 되겠습니다.
    # 예를 들어 이전엔 1 49 28... 이랬던 값들이 FALSE FALSE FALSE... 가 된다는 것이죠.
    # 이 말은 쉽게 말해서, bootstrap sample이 원 모집단의 j번째 값, 여기서는 4를
    # 포함하는지 아닌지 살펴보는 것입니다. 하나라도 TRUE가 있다면
    # bootstrap sample은 j 번째 값을 포함하고 있는 것이겠죠.
    
    # sum(sample(1:100, rep=TRUE)==4)
    # 이제 그 FALSE와 TRUE(TRUE가 있을지는 장담할 수 없습니다. 없을 수도 있겠죠.)로 이뤄진
    # 100개의 논리값들을 다 더합니다. 이 때 TRUE = 1, FALSE = 0 입니다.
    # 이 값이 0이면 bootstrap sample 안에는 j 번째 값이 없는 거고
    # 1 이상이면 한 번 이상 j 번째 값이 나왔다는 의미입니다.
    
    # sum(sample(1:100, rep=TRUE)==4) > 0
    # 다시 한 번 논리연산을 합니다. 위에서 논리값들을 다 더한 결과가 0보다 큰지 작은지를
    # 비교하여 그 결과를 store의 i 번째 공간에 저장한다는 뜻이지요.
    # store[i]는 TRUE 아니면 FALSE를 가집니다.
    # TRUE 는 bootstrap sample 에서 j 번째 값이 발견 되었다는 뜻이고
    # FALSE는 그렇지 않다는 뜻이겠죠.
}
# mean은 평균을 구하는 함수고, 아까 말한대로 TRUE = 1, FALSE = 0으로 간주하고 계산합니다.
# 아래 연산은 곧 bootstrap sampling을 만 번 했을 때, 그 중 j 번째 값이 하나 이상 들어있는
# sample들의 비율을 구하는 것이 되겠습니다. 결과가 대충 0.63 근처로 수렴할 것입니다.
mean(store)

0.6248

3. (a)

k-fold cross-validation이 어떻게 적용되는지 설명하자면,

1.Randomly divide the set of observations into k groups, or folds of approximately equal size.

2.The first fold is treated as a validation set, and the method is fit on the remaining k-1 folds.

3.The mean squared error, $MSE_1$, is then computed on the observations in the held-out fold.

4.This procedure is repeated k times.

5.Each time, a different group of observations is treated as a validation set.

6.This process results in k estimates of the test error.

7.The k-fold CV estimate is computed by averaging these values,

$$CV_{(k)} = \frac{1}{k}\sum_{i=1}^{k}\ MSE_i$$

자료를 크기가 서로 비슷한 k 개의 집단으로 나눈다. 첫 번째 집단을 validation set으로 삼고 나머지 집단들의 자료에 대해 fitting을 한다. fitting한 결과를 가지고 첫 번째 집단에 적용하여 MSE를 구한다. 이 과정을 마치면 validation set 을 두 번째 집단, 세 번째 집단... k 번째 집단으로 설정해가며 MSE를 각각 구한다. 이 MSE 들의 평균을 내는 것이 k-fold CV estimate 값이다.

3. (b)

k-fold CV가 validation set approach보다 좋은 점은

모든 data 들이 training set에 들어가기 때문에 overestimation의 위험을 줄일 수 있고

training set을 한 번만 구하는 게 아니라 여러번 구해서 mse의 평균을 내기 때문에 validation estimate의 분산을 더 줄일 수 있다.

LOOCV 보다 좋은 점은 calculation cost를 줄일 수 있다는 점이다.

4. 예측한 값의 표준편차를 추정하는 방법에 대해 자세히 설명할 것.

1번 문제에서 구한 $\alpha$의 추정치가 얼마나 정확한지 알 수 있을까? 추정치의 표준편차를 구하면 될 것이다.

그러나 그렇게 하려면 $\alpha$의 평균을 구하기 위해서 새로운 $X$, $Y$ 를 sampling을 해야 할 것이다. 이는 만만찮은 일이기 때문에, 새로이 $X$, $Y$를 구할 것이 아니라 이미 주어진 자료에서 bootstrap sample들을 구해서 $\alpha$의 추정치의 표준편차를 구하는 방법을 택한다.

원래의 자료로 주어진 집단을 $Z$ 라고 하고, bootstrap 방법을 통해 $Z$와 크기가 같은 sample 들을 취할 때, 이것을 각각 $Z^{*1}$, $Z^{*2}$, $Z^{*3}$,...,$Z^{*B}$ 라고 하자. B는 bootsrap 방법을 적용한 횟수이다.

그러면 $Z^{*1}$, $Z^{*2}$,...,$Z^{*B}$ 에 대하여 $\hat{\alpha}^{*1}$, $\hat{\alpha}^{*2}$,...,$\hat{\alpha}^{*B}$ 를 구할 수 있다.

이제 $\hat{\alpha}$의 표준오차를 구할 수 있다.

$$SE_{B}(\hat{\alpha}) = \sqrt{\frac{1}{B-1} \sum_{r=1}^{B} (\hat{\alpha}^{*r} - \frac{1}{B} \sum_{r'=1}^{B} \hat{\alpha}^{*r'})^2}$$

Applied

5. (a)

library(ISLR)

glm.fit <- glm(default~income+balance, data = Default, family='binomial')

summary(glm.fit)

Call:
glm(formula = default ~ income + balance, family = "binomial", 
    data = Default)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.4725  -0.1444  -0.0574  -0.0211   3.7245  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.154e+01  4.348e-01 -26.545  < 2e-16 ***
income       2.081e-05  4.985e-06   4.174 2.99e-05 ***
balance      5.647e-03  2.274e-04  24.836  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2920.6  on 9999  degrees of freedom
Residual deviance: 1579.0  on 9997  degrees of freedom
AIC: 1585

Number of Fisher Scoring iterations: 8

5. (b)

set.seed(1)

# 모든 과정을 통합해서 함수로 만들었습니다. 5(c) 에서 동일한 과정을 세 번 돌리라고 해서요.
five_b <- function(){
    # Split the sample set into a training set and a validation set
    # 이렇게 하면 정확하게 train과 test의 개수를 같게 할 수 있습니다.
    train<- sample(dim(Default)[1], dim(Default)[1]/2)

    # Fit a multiple logistic regression model using only the training observations
    glm.train <- glm(default~income+balance, Default, family='binomial',
                    subset = train)

    # Obtain a prediction of default status for each individual in the validation set
    # by computing the posterior probability of default for that individual, and
    # classifying the individual to the default category if the posterior probability
    # is greater than 0.5

    glm.probs <- predict(glm.train, Default[-train,], type = 'response')
    glm.pred <- rep("No", dim(Default)[1]/2)
    glm.pred[glm.probs>0.5] = "Yes"

    # Compute the validation set error, which is the fraction of the observations
    # in the validation set that are misclassified.
    return (mean(glm.pred != Default[-train,]$default))
}

five_b()

0.0286

5. (c)

five_b()

0.0236

five_b()

0.028

five_b()

0.0268

대략 validation set error가 0.026 정도에서 평균값을 가지는 것 같아요.

5. (d)

# 그냥 열 번 돌려봤습니다... 평균 내면 대충 0.026 근처의 값이 나옵니다.
for (i in 1:10){
    train<- sample(dim(Default)[1], dim(Default)[1]/2)

    glm.train <- glm(default~income+balance+student, Default, family='binomial',
                    subset = train)

    glm.probs <- predict(glm.train, Default[-train,], type = 'response')
    glm.pred <- rep("No", dim(Default)[1]/2)
    glm.pred[glm.probs>0.5] = "Yes"

    print(mean(glm.pred != Default[-train,]$default))
    
}

[1] 0.0276
[1] 0.0258
[1] 0.026
[1] 0.0268
[1] 0.0264
[1] 0.0258
[1] 0.0258
[1] 0.0282
[1] 0.0252
[1] 0.026

student를 dummy variable로 집어 넣어도 validation set error에는 별 영향이 없어 보여요.

6. (a)

# Using the summary()...
glm.train <- glm(default~income+balance, Default, family='binomial')

# ... and glm() functions, determine the estimated standard errors for the
# coefficients associated with income and balance in 
# a multiple logistic regression model that uses both predictors.
summary(glm.train)

Call:
glm(formula = default ~ income + balance, family = "binomial", 
    data = Default)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.4725  -0.1444  -0.0574  -0.0211   3.7245  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.154e+01  4.348e-01 -26.545  < 2e-16 ***
income       2.081e-05  4.985e-06   4.174 2.99e-05 ***
balance      5.647e-03  2.274e-04  24.836  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2920.6  on 9999  degrees of freedom
Residual deviance: 1579.0  on 9997  degrees of freedom
AIC: 1585

Number of Fisher Scoring iterations: 8

income의 경우에는 표준오차가 0.000004985, balance의 경우에는 0.0002274가 되겠습니다.

6. (b)

boot.fn <- function(data, index){
    
    model <- glm(default~income+balance, data = data, family='binomial', 
                 subset=index)
    return (coef(model))
    
}

6. (c)

library(boot)

boot(Default, boot.fn, 1000)

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Default, statistic = boot.fn, R = 1000)


Bootstrap Statistics :
         original        bias     std. error
t1* -1.154047e+01 -7.844862e-03 4.245605e-01
t2*  2.080898e-05  5.396354e-08 4.587088e-06
t3*  5.647103e-03  2.234814e-06 2.269421e-04

위의 결과표는 위에서부터 각각 절편, income, balance에 해당하는 coefficient들과 그 값들의 표준오차를 나타냅니다.

6. (d)

glm을 사용한 결과와 bootstrapping을 한 결과를 놓고 비교해보면,

glm 사용 시 계수들의 표준오차는 4.348e-01, 4.985e-06, 2.274e-04 이고, boot 사용 시에는 4.245e-01, 4.587e-06, 2.269e-04 입니다.

결론은 뭐... 비슷해요.

7. (a)

glm.fit <- glm(Direction~Lag1+Lag2, data=Weekly, family='binomial')
summary(glm.fit)

Call:
glm(formula = Direction ~ Lag1 + Lag2, family = "binomial", data = Weekly)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.623  -1.261   1.001   1.083   1.506  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.22122    0.06147   3.599 0.000319 ***
Lag1        -0.03872    0.02622  -1.477 0.139672    
Lag2         0.06025    0.02655   2.270 0.023232 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1496.2  on 1088  degrees of freedom
Residual deviance: 1488.2  on 1086  degrees of freedom
AIC: 1494.2

Number of Fisher Scoring iterations: 4

7. (b)

sub <- 2:nrow(Weekly)
glm.fit2 <- glm(Direction~Lag1+Lag2, data=Weekly, family='binomial', subset=sub)
summary(glm.fit2)

Call:
glm(formula = Direction ~ Lag1 + Lag2, family = "binomial", data = Weekly, 
    subset = sub)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6258  -1.2617   0.9999   1.0819   1.5071  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.22324    0.06150   3.630 0.000283 ***
Lag1        -0.03843    0.02622  -1.466 0.142683    
Lag2         0.06085    0.02656   2.291 0.021971 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1494.6  on 1087  degrees of freedom
Residual deviance: 1486.5  on 1085  degrees of freedom
AIC: 1492.5

Number of Fisher Scoring iterations: 4

7. (c)

predict(glm.fit2, Weekly[1,], type='response')

1: 0.571392320520443

문제에 주어진 조건에 따라 $P(Direction= Up$ | $Lag1, Lag2) > 0.5$ 이면, 위의 결과 값은 첫 번째 관찰값에 대해 "Up" 이라고 예측한 것입니다 (P = 0.571392). 그러니 예측을 실패한 것이 되겠습니다.

7. (d)

error_count <- rep(NA, nrow(Weekly))
for (i in 1:nrow(Weekly)){
    
    data = Weekly[-i, ]
    glm.fit_i <- glm(Direction~Lag1+Lag2, data=data, family='binomial')
    glm.prob <- predict(glm.fit_i, Weekly[i,], type='response')
    if(glm.prob > 0.5)
        pred <- "Up"
    else
        pred <- "Down"
    if (pred != Weekly$Direction[i])
        error_count[i] <- 1
    else
        error_count[i] <- 0
}

0.449954086317723

7. (e)

mean(error_count)

0.449954086317723

LOOCV의 error rate는 0.45가 되겠습니다. 생각보다 많이 못 맞추네요.

8. (a)

set.seed(1)
y = rnorm(100)
x = rnorm(100)
y = x-2*x^2+rnorm(100)

위 data set 에서 n은 100이고, p는 2가 되겠습니다. 이걸 방정식으로 써 보면

$$Y = X-2X^2 + \epsilon$$

와 같이 쓸 수 있겠습니다.

8. (b)

plot(x, y)

X와 Y의 관계가 linear 하지 않고 quadratic 한 것을 알 수 있습니다.

8. (c)

dataset = data.frame(x, y)

# i의 숫자에 따라 seed를 다르게 정하는 방법을 썼습니다.
eight_c <- function(i){
    
    set.seed(i)
    
    # i.
    model_1 <- glm(y~x, data = dataset)
    print(cv.glm(dataset, model_1)$delta)

    # ii.
    model_2 <- glm(y~x+I(x^2), data = dataset)
    print(cv.glm(dataset, model_2)$delta)

    # iii.
    model_3 <- glm(y~x+I(x^2)+I(x^3), data = dataset)
    print(cv.glm(dataset, model_3)$delta)

    # iv.
    model_4 <- glm(y~x+I(x^2)+I(x^3)+I(x^4), data = dataset)
    print(cv.glm(dataset, model_4)$delta)
    
}

eight_c(1)

[1] 5.890979 5.888812
[1] 1.086596 1.086326
[1] 1.102585 1.102227
[1] 1.114772 1.114334

8. (d)

eight_c(5)

[1] 5.890979 5.888812
[1] 1.086596 1.086326
[1] 1.102585 1.102227
[1] 1.114772 1.114334

다른 random seed를 정해도 값이 똑같아요. 어차피 cv.glm 함수는 주어진 자료를 K 개로 나눠서 살펴보는데, K가 주어지지 않을 때에는 K = n, 곧 observation의 수와 같기 때문에 n 개의 자료를 n 개로 나누면 single observation들을 매번 살펴보는 것과 같기 때문입니다.

즉, seed를 어떻게 정하든 상관 없이, 각 test set에 들어가는 값은 항상 하나의 observation이고 training set에 들어가는 값은 나머지 observations이기 때문에 seed를 바꾸는 것은 결과 값에 영향을 끼치지 않을 것이라는 이야기입니다.

8. (e)

두 번째 모델이 가장 error rate가 적습니다. 당연히 원래 x와 y는 quadratic한 관계를 지니기 때문에 2차 항을 넣은 모델이 가장 error rate가 적겠지요.

8. (f)

summary(glm(y~x))

Call:
glm(formula = y ~ x)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-7.3469  -0.9275   0.8028   1.5608   4.3974  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -1.8185     0.2364  -7.692 1.14e-11 ***
x             0.2430     0.2479   0.981    0.329    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 5.580018)

    Null deviance: 552.21  on 99  degrees of freedom
Residual deviance: 546.84  on 98  degrees of freedom
AIC: 459.69

Number of Fisher Scoring iterations: 2

summary(glm(y~x+I(x^2)))

Call:
glm(formula = y ~ x + I(x^2))

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.89884  -0.53765   0.04135   0.61490   2.73607  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.09544    0.13345  -0.715    0.476    
x            0.89961    0.11300   7.961 3.24e-12 ***
I(x^2)      -1.86665    0.09151 -20.399  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 1.06575)

    Null deviance: 552.21  on 99  degrees of freedom
Residual deviance: 103.38  on 97  degrees of freedom
AIC: 295.11

Number of Fisher Scoring iterations: 2

summary(glm(y~x+I(x^2)+I(x^3)))

Call:
glm(formula = y ~ x + I(x^2) + I(x^3))

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.87250  -0.53881   0.02862   0.59383   2.74350  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.09865    0.13453  -0.733    0.465    
x            0.95551    0.22150   4.314  3.9e-05 ***
I(x^2)      -1.85303    0.10296 -17.998  < 2e-16 ***
I(x^3)      -0.02479    0.08435  -0.294    0.769    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 1.075883)

    Null deviance: 552.21  on 99  degrees of freedom
Residual deviance: 103.28  on 96  degrees of freedom
AIC: 297.02

Number of Fisher Scoring iterations: 2

summary(glm(y~x+I(x^2)+I(x^3)+I(x^4)))

Call:
glm(formula = y ~ x + I(x^2) + I(x^3) + I(x^4))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.8914  -0.5244   0.0749   0.5932   2.7796  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.13897    0.15973  -0.870 0.386455    
x            0.90980    0.24249   3.752 0.000302 ***
I(x^2)      -1.72802    0.28379  -6.089  2.4e-08 ***
I(x^3)       0.00715    0.10832   0.066 0.947510    
I(x^4)      -0.03807    0.08049  -0.473 0.637291    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 1.084654)

    Null deviance: 552.21  on 99  degrees of freedom
Residual deviance: 103.04  on 95  degrees of freedom
AIC: 298.78

Number of Fisher Scoring iterations: 2

summary(glm(y~poly(x,4)))

Call:
glm(formula = y ~ poly(x, 4))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.8914  -0.5244   0.0749   0.5932   2.7796  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -1.8277     0.1041 -17.549   <2e-16 ***
poly(x, 4)1   2.3164     1.0415   2.224   0.0285 *  
poly(x, 4)2 -21.0586     1.0415 -20.220   <2e-16 ***
poly(x, 4)3  -0.3048     1.0415  -0.293   0.7704    
poly(x, 4)4  -0.4926     1.0415  -0.473   0.6373    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 1.084654)

    Null deviance: 552.21  on 99  degrees of freedom
Residual deviance: 103.04  on 95  degrees of freedom
AIC: 298.78

Number of Fisher Scoring iterations: 2

각 모델을 I를 사용해서 보거나, 아니면 poly를 써서 보더라도 1차항과 2차항은 통계적으로 유의미하고, 나머지 항들은 통계적으로 유의미하지 않습니다. 이는 cross-validation 결과와 일치합니다.

9. (a)

library(MASS)
summary(Boston)

      crim                zn             indus            chas        
 Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46   Min.   :0.00000  
 1st Qu.: 0.08204   1st Qu.:  0.00   1st Qu.: 5.19   1st Qu.:0.00000  
 Median : 0.25651   Median :  0.00   Median : 9.69   Median :0.00000  
 Mean   : 3.61352   Mean   : 11.36   Mean   :11.14   Mean   :0.06917  
 3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10   3rd Qu.:0.00000  
 Max.   :88.97620   Max.   :100.00   Max.   :27.74   Max.   :1.00000  
      nox               rm             age              dis        
 Min.   :0.3850   Min.   :3.561   Min.   :  2.90   Min.   : 1.130  
 1st Qu.:0.4490   1st Qu.:5.886   1st Qu.: 45.02   1st Qu.: 2.100  
 Median :0.5380   Median :6.208   Median : 77.50   Median : 3.207  
 Mean   :0.5547   Mean   :6.285   Mean   : 68.57   Mean   : 3.795  
 3rd Qu.:0.6240   3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.: 5.188  
 Max.   :0.8710   Max.   :8.780   Max.   :100.00   Max.   :12.127  
      rad              tax           ptratio          black       
 Min.   : 1.000   Min.   :187.0   Min.   :12.60   Min.   :  0.32  
 1st Qu.: 4.000   1st Qu.:279.0   1st Qu.:17.40   1st Qu.:375.38  
 Median : 5.000   Median :330.0   Median :19.05   Median :391.44  
 Mean   : 9.549   Mean   :408.2   Mean   :18.46   Mean   :356.67  
 3rd Qu.:24.000   3rd Qu.:666.0   3rd Qu.:20.20   3rd Qu.:396.23  
 Max.   :24.000   Max.   :711.0   Max.   :22.00   Max.   :396.90  
     lstat            medv      
 Min.   : 1.73   Min.   : 5.00  
 1st Qu.: 6.95   1st Qu.:17.02  
 Median :11.36   Median :21.20  
 Mean   :12.65   Mean   :22.53  
 3rd Qu.:16.95   3rd Qu.:25.00  
 Max.   :37.97   Max.   :50.00  

mean(Boston$medv)

22.5328063241107

9. (b)

sqrt(var(Boston$medv)/length(Boston$medv))

0.408861147497535

9. (c)

fn <- function(data, index){
    
    return (mean(data[index]))
}

boot(Boston$medv, fn, 1000)

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Boston$medv, statistic = sample_error, R = 1000)


Bootstrap Statistics :
    original      bias    std. error
t1* 22.53281 0.003793281    0.403674

결과 값은 b에서 구한 것과 비슷합니다(0.40886.. vs 0.403674)

9. (d)

ci <- c(22.53281-1.96*0.403674, 22.53281+1.96*0.403674)
ci

1. 21.74160896
2. 23.32401104

t.test(Boston$medv)

	One Sample t-test

data:  Boston$medv
t = 55.111, df = 505, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 21.72953 23.33608
sample estimates:
mean of x 
 22.53281 

bootstrap 결과를 통해 구한 95% 신뢰구간은 t.test를 통해서 구한 신뢰구간과 그렇게 큰 차이를 보이지 않습니다.

[21.74160896, 23.32401104] vs [21.72953, 23.33608]

9. (e)

median(Boston$medv)

21.2

median의 표본평균 추정치는 21.2 입니다.

9. (f)

fn2 <- function(data, index){
    return (median(data[index]))
}

boot(Boston$medv, fn2, 1000)

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Boston$medv, statistic = fn2, R = 1000)


Bootstrap Statistics :
    original   bias    std. error
t1*     21.2 -0.00475   0.3827354

median의 표본평균 추정치 21.2에 대하여 표준오차는 0.3827354 입니다.

9. (g)

quantile(Boston$medv, c(0.1))

10%: 12.75

9. (f)

fn3 <- function(data, index){
    return (quantile(data[index], c(0.1)))
}

boot(Boston$medv, fn3, 1000)

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Boston$medv, statistic = fn3, R = 1000)


Bootstrap Statistics :
    original   bias    std. error
t1*    12.75 -0.00435    0.511498

$\hat\mu_{0.1}$ 의 표준오차는 0.511498이 되겠습니다.